function D_matrixs(E::Float64, μ::Float64)
	D0 = E/(1.0 - μ*μ)
	return D0*[1.0 μ 0.0; μ 1.0 0.0; 0.0 0.0 0.5*(1.0 - μ)]
end

function D_matrixc(E::Float64, G::Float64)
	return [E 0.0; 0.0 5.0*G/6.0]
end

function ksteel(xi::Float64, xj::Float64)
   k = zeros(Float64, 30, 30)
   det = (xj - xi)/2.0

   a1 = -0.5*hweb - fbot; b1 = -0.5*hweb;
   a2 = -0.5*hweb; b2 = 0.5*hweb;
   a3 = 0.5*hweb; b3 = 0.5*hweb + ftop;
   D0 = D_matrixs(Es, μs)

   # I-shaped steel beam integration
   # part 1, lower flange
   for i = 1:ngpx # integrate in x axis
        for j = 1:ngpz # integrate in z axis
            ξ = gaussdatax[i, 1]
            z = 0.5*(b1 - a1)*gaussdataz[j, 1] + 0.5*(b1 + a1)
            w = gaussdatax[i, 2]*0.5*(b1 - a1)*gaussdataz[j, 2]
			B = B_matrix_s(ξ, z, xi, xj)
			int = (B') * D0 * B
            k = k + int*det*wbot*w
        end
    end
	# part 2, web
	for i = 1:ngpx # integrate in x axis
		 for j = 1:ngpz # integrate in z axis
			 ξ = gaussdatax[i, 1]
			 z = 0.5*(b2 - a2)*gaussdataz[j, 1] + 0.5*(b2 + a2)
			 w = gaussdatax[i, 2]*0.5*(b2 - a2)*gaussdataz[j, 2]
			 B = B_matrix_s(ξ, z, xi, xj)
			 int = (B') * D0 * B
			 k = k + int*det*wweb*w
		 end
	 end
	 # part 3, upper flange
	 for i = 1:ngpx # integrate in x axis
		  for j = 1:ngpz # integrate in z axis
			  ξ = gaussdatax[i, 1]
			  z = 0.5*(b3 - a3)*gaussdataz[j, 1] + 0.5*(b3 + a3)
			  w = gaussdatax[i, 2]*0.5*(b3 - a3)*gaussdataz[j, 2]
			  B = B_matrix_s(ξ, z, xi, xj)
			  int = (B') * D0 * B
			  k = k + int*det*wtop*w
		  end
	  end
	return k
end

function kcs(xi::Float64, xj::Float64)
	k = zeros(30,30)
	det = (xj - xi)/2.0

	for i = 1:ngpx
			ξ = gaussdatax[i, 1]
			w = gaussdatax[i, 2]
			n1 = N_ucs(ξ, xi, xj, h1, h2)
			n2 = N_wcs(ξ, xi, xj, h1, h2)
			int = n1' * kucs * n1 + n2' * kwcs * n2
			k = k + int*det*w
	end
	return k
end

function krebar(xi::Float64, xj::Float64, z::Float64)
	# z: rebar possition in o-zs coordinate
	k = zeros(30, 30)
	ngpx = size(gaussdatax,1)
	det = (xj - xi)/2.0
	for i = 1:ngpx
		ξ = gaussdatax[i, 1]
		w = gaussdatax[i, 2]
		n = N_ϵcx(ξ, z, xi, xj)
		int = n' * n
		k = k + int*Ar*(Er-Ec)*det*w
	end
	return k
end

function kconcrete(xi::Float64, xj::Float64)
	k = zeros(30, 30)
	det = (xj - xi)/2.0
	D0 = D_matrixc(Ec, Gc)
	for i = 1:ngpx
		for j = 1:ngpz
			ξ = gaussdatax[i, 1]
			z = 0.5*Hc*gaussdataz[j, 1]
			w = gaussdatax[i, 2]*0.5*Hc*gaussdataz[j, 2]
			B = B_matrix_c(ξ, z, xi, xj)
			int = (B') * D0 * B
			k = k + int*det*Bc*w
		end
	end
	return k
end

function R0(xi::Float64, xj::Float64, q::Float64)
	r = zeros(ndof)
	L = xj - xi
	r[[3,13,23]] = [(L*q)/6, (2*L*q)/3, (L*q)/6]
	return r
end

function R1(nel, nt)
	int = zeros(ndof);
	num = gnum[nel,:]
	elemcoord = gcoord[num,:]
	xi = elemcoord[1]; xj = elemcoord[3]
	det = (xj - xi)/2.0
	for i = 1:ngpz
		for j = 1:ngpx
			z = 0.5*Hc*gaussdataz[i, 1]
			ξ = gaussdatax[j, 1]
			w = gaussdatax[j, 2]*0.5*Hc*gaussdataz[i, 2]
			s = 0.0
			for i = 1:(nt-1)
				s = s + psi(nt, i, time0)*dot(StressHistory[i, nel, :], Nstress(ξ, z))
			end
			Ks = reshape(N_ϵcx(ξ, z, xi, xj), ndof)*(en(nt, time0)*1e6*epsilon_shn(nt, time0) - s)
			int = int + Bc*det*w*Ks
		end
	end
	for j = 1:ngpx
		ξ = gaussdatax[j, 1]
		w = gaussdatax[j, 2]
		s = 0.0
		for i = 1:(nt-1)
			s = s + psi(nt, i, time0)*dot(StressHistory[i, nel, :], Nstress(ξ, zr))
		end
		Ks = reshape(N_ϵcx(ξ, zr, xi, xj), ndof)*(s - en(nt, time0)*1e6*epsilon_shn(nt, time0))
		int = int + Ar*det*w*Ks
	end

	for i = 1:ngpz
		for j = 1:ngpx
			z = 0.5*Hc*gaussdataz[i, 1]
			ξ = gaussdatax[j, 1]
			w = gaussdatax[j, 2]*0.5*Hc*gaussdataz[i, 2]
			s = 0.0
			for i = 1:(nt-1)
				s = s + psi(nt, i, time0)*dot(StressHistorytao[i, nel, :], Nstress(ξ, z))
			end
			Ks = -reshape(N_γc(ξ, xi, xj), ndof)*s
			int = int + Bc*det*w*Ks
		end
	end

	return int
end

function Nstress(ξ, z)
	nx = lagrange(gx, ξ)
	nz = lagrange(0.5*Hc*gz, z)
	np = zeros(9)
	count = 0
	for iz = 1:3
		for ix = 1:3
			count = count + 1
			np[count] = nx[ix]*nz[iz]
		end
	end
	return np
end

function strainc(nel, sol, ξ, z)
   elemcoord = gcoord[gnum[nel,:],:]
   xi = elemcoord[1]
   xj = elemcoord[3]
   de = [extract(sol, g_g[k, nel]) for k = 1:ndof]
   return dot(N_ϵcx(ξ, z, xi, xj), de)
end

function strainctao(nel, sol, ξ, z)
   elemcoord = gcoord[gnum[nel,:],:]
   xi = elemcoord[1]
   xj = elemcoord[3]
   de = [extract(sol, g_g[k, nel]) for k = 1:ndof]
   return dot(N_γc(ξ, xi, xj), de)
end

function stressc(sol, nel, ξ, z, nt)
	#(* sol is the solution of the ntth time tick *)
	s = 0.0
	for i = 1:(nt-1)
		s = s + psi(nt, i, time0)*dot(StressHistory[i, nel,:], Nstress(ξ, z))
	end
	r = en(nt, time0)*1e6*(strainc(nel, sol, ξ, z) - epsilon_shn(nt, time0)) + s
	return r
end

function stressctao(sol, nel, ξ, z, nt)
	#(* sol is the solution of the ntth time tick *)
	s = 0.0
	for i = 1:(nt-1)
		s = s + psi(nt, i, time0)*dot(StressHistorytao[i, nel,:], Nstress(ξ, z))
	end
	gc = (5.0/6.0) * en(nt, time0)/(2.0*(1 + μc)) * 1e6
	r = gc * strainctao(nel, sol, ξ, z) + s
	return r
end

function record(sol, nt)
	#   (* sol is the solution of the ntth time tic *)
	sel = zeros(nels, 9)
	seltao = zeros(nels, 9)
	for nel = 1:nels
		count = 0
		for iz = 1:3
			for ix = 1:3
				count = count + 1
				z = 0.5*Hc*gz[iz]
				ξ = gx[ix]
				sel[nel, count] = stressc(sol, nel, ξ, z, nt)
				seltao[nel, count] = stressctao(sol, nel, ξ, z, nt)
			end
		end
	end
	StressHistory[nt,:,:] = sel
	StressHistorytao[nt,:,:] = seltao
	return nothing
end

function assemble0(gnum::Array{Int64,2}, gcoord::Array{Float64,2}, nf::Array{Int64,2}, nels::Int64, neq::Int64)
    gK = zeros(Float64, neq, neq)
    gR = zeros(Float64, neq)

    for iel = 1:nels
        num = gnum[iel,:]
        elemcoord = gcoord[num,:]
		xi = elemcoord[1]; xj=elemcoord[3];
        g = numtog(num, nf)
        kma = ksteel(xi, xj) + kcs(xi,xj)
		rvec = R0(xi, xj, q0)
	    for i = 1:length(g)
			row = g[i]
			if row != 0
				gR[row] = gR[row] + rvec[i]
			end
	        for j = 1:length(g)
	            col = g[j]
	            if row*col != 0
	                gK[row, col] = gK[row, col] + kma[i, j]
	            end
	        end
	    end
	end

	return (gK, gR)
end

function assemble1(gnum::Array{Int64,2}, gcoord::Array{Float64,2}, nf::Array{Int64,2}, nels::Int64, neq::Int64)
    gK = zeros(Float64, neq, neq)
    gR = zeros(Float64, neq)

    for iel = 1:nels
        num = gnum[iel,:]
        elemcoord = gcoord[num,:]
		xi = elemcoord[1]; xj=elemcoord[3];
        g = numtog(num, nf)
        kma = krebar(xi,xj,zr) + kconcrete(xi,xj)
		rvec = R1(iel, nt)
	    for i = 1:length(g)
			row = g[i]
			if row != 0
				gR[row] = gR[row] + rvec[i]
			end
	        for j = 1:length(g)
	            col = g[j]
	            if row*col != 0
	                gK[row, col] = gK[row, col] + kma[i, j]
	            end
	        end
	    end
	end

	return (gK, gR)
end

function assemble_loads0(gnum::Array{Int64,2}, gcoord::Array{Float64,2}, nf::Array{Int64,2}, nels::Int64, neq::Int64)
    gR = zeros(Float64, neq)
    for iel = 1:nels
        num = gnum[iel,:]
        elemcoord = gcoord[num,:]
		xi = elemcoord[1]; xj=elemcoord[3]
        g = numtog(num, nf)
		rvec = R0(xi, xj, q0)
	    for i = 1:length(g)
			row = g[i]
			if row != 0
				gR[row] = gR[row] + rvec[i]
			end
	    end
	end

	return gR
end

function assemble_loads1(gnum::Array{Int64,2}, gcoord::Array{Float64,2}, nf::Array{Int64,2}, nels::Int64, neq::Int64)
    gR = zeros(Float64, neq)
    for iel = 1:nels
        num = gnum[iel,:]
        g = numtog(num, nf)
		rvec = R1(iel, nt)
	    for i = 1:length(g)
			row = g[i]
			if row != 0
				gR[row] = gR[row] + rvec[i]
			end
	    end
	end
	return gR
end
